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We have studied the three-dimensional Ising spin glass with a ±J distribution by Monte Carlo 
simulations. Using larger sizes and much better statistics than in earlier work, a finite size scaling 
analysis shows quite strong evidence for a finite transition temperature, T c , with ordering below T c . 
Our estimate of the transition temperature is rather lower than in earlier work, and the value of the 
correlation length exponent, v, is somewhat higher. Because there may be (unknown) corrections 
to finite size scaling, we do not completely rule out the possibility that T c = or that T c is finite 
but with no order below T c . However, from our data, these possibilities seem less likely. 
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The question of whether there is a finite transition tem- 
perature, T c , in an Ising spin glass in three dimensions 
has aroused a lot of interest for the last two decades 1 , 
and the consensus of opinion has changed several times. 
About one decade ago, several pieces of work 2-5 seemed 
to show that there is a finite T c , and this conclusion has 
generally been restated since then 6 . However, on closer 
inspection, the case is not completely closed. For exam- 
ple, the work of one of us 2 , henceforth referred to as BY, 
is unable to rule out the possibility that T c = and the 
correlation length, £, diverges exponentially as T — > 0, as 
happens in the two-dimensional Heisenberg ferromagnet. 
The data is also consistent with a line of critical points 
terminating at T c ~ 1.2, as occurs in the Kosterlitz- 
Thouless-Berezinskii theory of the two-dimensional XY 
ferromagnet. In this scenario there would be no long 
range spin glass order below T c . Furthermore, recent 
results of Marinari et al. 7 were found to be consistent 
both with a finite T c and with a zero temperature transi- 
tion where the correlation length diverges exponentially, 
£ ~ exp(A/T 4 ). We therefore feel there are three possible 
scenarios, consistent with existing work: 

(i) T c is finite and there is spin glass order at lower 
temperatures, 



ii 



T c is finite but there is a line of critical points (i.e. 
no spin glass order) at lower temperatures, 



(iii) T c = and the correlation length diverges expo- 
nentially as T — > 0. 

During the last decade available computer power has in- 
creased enormously so, given these uncertainties, it is 
useful to look at the problem again. The calculations 
presented here are similar to those of BY, but we are able 
to study larger system sizes in the temperature range of 
interest and obtain much better statistics by averaging 
over many more samples. As a result, unlike BY, we are 
able to see clear evidence for ordering below a finite T c . 
The Hamiltonian is 



TABLE I. For each size, L, we show the largest value of to, 
(where, as explained in the text, the simulation ran for 3 to 
sweeps) and the minimum number of samples, N s . 
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largest to 


minimum N s 
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4 x 10 5 


8192 


8 


1 x 10 6 


8192 


12 


8 x 10 6 


6880 


16 


15 x 10 6 


3392 


24 


5 x 10 6 


2080 



ft — JijSiSj 



(1) 



where the spins Si take values ±1, and the nearest neigh- 
bor interactions, Jij take values ±1 with equal proba- 
bility. The simple cubic lattice contains N = L 3 spins 
and has periodic boundary conditions. In some previ- 
ous work, {Jij} was generated so that the the number 
of ferromagnetic couplings is exactly the same as that of 
antiferromagnetic couplings. We do not impose such a 
condition in the present work. 

The Monte Carlo simulation uses a multispin coding 
technique 8 in which each spin and bond is represented 
by a single bit of a computer word. On a 32 bit machine 
we then flip in parallel 32 spins (on the same lattice site 
but in different samples with different realizations of the 
disorder). For this method to be efficient the same ran- 
dom number is used for each bit 9 . We use a shift reg- 
ister random number generator 10,11 , commonly known 
as R250. The code runs at 27 million spin updates per 
second on one node (IBM 390 RISC workstation) of the 
SP2 computer at the Maui High Performance Comput- 
ing Center. Since we need many more than 32 samples, 
we ran the same code independently on many nodes at 
the same time. Each node produces its own output file 
from which the final averaging is easily done using a unix 
shell script. Monte Carlo simulations of random systems 
thus provide an example where parallel computing can 
be done in a trivial (and almost perfectly efficient!) way. 
The total CPU time used for the data presented here is 
about 9 node-years. 
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To get good statistics we average over a large number 
of samples, N s , where for each size, N s is at least the 
value in the third column in table I. After to sweeps for 
equilibration, an additional 2to sweeps are carried out for 
measurements. For each size, the largest value of to used 
is also shown in table I (this is for the lowest temperature: 
at higher temperatures many fewer sweeps are generally 
needed). 

As usual 2 , for each realization of the bonds, two copies 
of the system are studied with different initial values of 
the spins and different random numbers for generating 
the spin flips. Of particular importance is the overlap 
between the two copies, 

9 = ^XX^ 2) > ( 2 ) 

8 = 1 

where the labels "1" and "2" denote the copies. From 
measurements of q we compute the Binder ratio 12,2 



where the average (...) denotes both a thermal average for 
a given set of bonds and an average over the disorder 13 . 
At high temperature, g — > 0, whereas g — » 1 in the spin 
glass phase, at least if there is a unique thermodynamic 
state. 
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FIG. 1. Results for the Binder ratio g, defined in Eq. (3), 
for different sizes and temperatures. The lines are smooth 
curves through the data and are only intended as guides to 
the eye. 

Because g is dimensionless it has the finite size scaling 
form 2 

g = g(L^(T-T c )) (4) 

and so is independent of L at T c . The behavior of g is 
different for each of the three scenarios discussed above: 



(i) the curves for g will intersect at T c and splay out 
again at lower T (with the larger sizes having the 
larger values, the opposite of the situation above 
Tc), 

(ii) the curves for g will come together at T c and then 
stick together at lower T, 

(iii) the curves will merge once £ ^> L, but data for 
larger sizes will merge to this common curve at suc- 
cessively lower temperatures. 
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FIG. 2. An enlarged view of the data in Fig. 1 in the crucial 
region where the curves come together. 

In addition to g, we also computed the spin glass sus- 
ceptibility, 

Xsa = N{q 2 ) , (5) 

and P{q), the distribution of q. These have the finite size 
scaling forms, 2 

XsG (^(T-T c )) , (6) 

and 

P(q) = L^P{L^q, I}> V {T - T c )) , (7) 

where [3 is the order parameter exponent and is related 
to rj, which gives the power law decay of the correlations 
at the critical point, by 

£ = l(d-2 + T,). (8) 

It is very important to ensure that enough Monte Carlo 
sweeps have been carried out to equilibrate the sample. 
Following BY we compare the results for g obtained, as 
described above, from the overlap between two replicas 
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with the results obtained from one replica at two different 
times (see BY for details). BY found that these two 
estimates approach the equilibrium value from opposite 
directions as the length of the simulation increased. Once 
the two values agree, they do not change further if more 
sweeps are carried out. We have also tested this by doing 
the run for L = 8, T = 0.9618 for an order of magnitude 
longer time than needed for the two estimates to agree. 
Again we find that there is no subsequent change within 
our (much smaller) errors. 
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FIG. 3. A scaling plot for g according to the form in Eq. (4). 

Our data for g is shown in Fig. 1 and an enlarged 
view of the region where the curves for different sizes 
intersect is shown in Fig. 2. From Fig. 2 one sees clear 
evidence for splaying out of the data below a temperature 
of about 1.10. Estimating T c to be approximately 1.10 
from the intersection point we can scale most of the data 
according to Eq. (4) with v = 2.0, see Fig. 3. The only 
point which does not lie on a common curve is the result 
for L = 24, T = 1.1948, which is significantly higher. 
One can see from Fig. 1 that this point has almost the 
same value of g as the data for L = 16 at the same 
temperature. This data point being rather higher than 
expected may reflect corrections to finite size scaling, and 
indicate that the true critical temperature is higher than 
the straightforward estimate based on data for g with 
L < 16. 

Once T c has been estimated, one can obtain (3/v, or 
equivalently rj, from the expected scaling form of P(q) at 
criticality given by Eq. (7) with T = T c . The correspond- 
ing plot is shown in Fig. 4 for T = 1.1113 (well within 
the bounds of our estimate of T c ), and has f3/v = 0.3 
which corresponds to rj = —0.4 from Eq. (8) with d = 3. 

We have also performed finite size scaling plots for Xsa 
according to Eq. (6). This data does not locate T c pre- 
cisely, so we have used the same T c as obtained from the 
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FIG. 4. A scaling plot for P(q) at T = 1.1113 (which is 
close to the the critical point) according to the form in Eq. (7). 
According to Eq. (8), the value jijv = 0.3 corresponds to 
■q = -0.4. 

scaling plot for g in Fig. 3, i.e. T c = 1.10. Furthermore 
the value of rj is constrained by requiring that the data 
scales at T c and from Fig. 4 this gives rj = —0.4 The 
only remaining parameter is v and the best fit, shown in 
Fig. 5, is for v = 1.6. 

The values for v obtained from g and Xsa are some- 
what different. If we try to use v = 2.0 in the data for 
Xsa or v = 1.6 in the data for g, the fit is visibly worse. 
Presumably this difference indicates that corrections to 
finite size scaling are not negligible for the range of sizes 
that we can study. Taking into account all the data we 
estimate 

T c = 1.11 ±0.04 
v = 1.7 ± 0.3 

r, = -0.35 ± 0.05 . (9) 

As discussed above, the L = 24 data indicates that T c 
may be higher than that estimated from the intersections 
of g for L < 16. This is reflected in the estimated error 
for T c in Eq. (9). The estimated errors in v and rj then 
come largely from the uncertainty in T c . Our value of 
T c is rather lower than earlier estimates which were close 
to 1.2, and the value of v is higher, previous estimates 
generally being in the vicinity of 1.3. Our value of rj is 
not very different from earlier estimates. 

To conclude, we have found evidence for a finite transi- 
tion temperature with spin glass order below T c , scenario 
(i) above. However, it is difficult to estimate the size 
of systematic errors, such as possible correlations in the 
random numbers (though we believe that these are very 
small 11 ), and corrections to finite size scaling. Because 
of this, and because the crossing of the data for g that 
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FIG. 5. A scaling plot for %sg according to the form in 
Eq. (6). 

we observe in Fig 2 is rather small, we cannot rule out, 
for sure, the other two possibilities, i.e. scenario (ii) in 
which T c is finite but there is no spin glass order at lower 
temperature, or scenario (iii) in which T c = 0. However, 
from our data, these possibilities now seem less likely. 

Since the present study required a substantial com- 
puter effort, an investigation of larger sizes, which is nec- 
essary to confirm scenario (i) beyond reasonable doubt, 
may need a better algorithm than single spin flip Monte 
Carlo. There are already some promising results from 
the "replica exchange" method 14 (where, in addition to 
local moves, global moves are made which cause the tem- 
perature of the system to cycle up and down). 
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